home *** CD-ROM | disk | FTP | other *** search
/ Visual Basic Source Code / Visual Basic Source Code.iso / vbsource / optivc32 / mdstd.h < prev    next >
Encoding:
C/C++ Source or Header  |  1999-03-06  |  57.4 KB  |  1,086 lines
  1. /*  MDstd.h
  2.  
  3.   matrix management functions:
  4.   manipulations on matrices of data type "double"
  5.   (double-precision real numbers)
  6.  
  7.   Copyright (c) 1996-1999 by Martin Sander
  8.   All Rights Reserved.
  9. */
  10.  
  11. #if !defined( __MDSTD_H )
  12. #define __MDSTD_H
  13. #if !defined( __MATLIB_H )
  14. #include <MatLib.h>
  15. #endif
  16. #if !defined( __VDSTD_H )
  17. #include <VDstd.h>
  18. #endif
  19. #if !defined( __VDMATH_H )
  20. #include <VDmath.h>
  21. #endif
  22.  
  23. #ifdef __cplusplus
  24. extern "C" {
  25. #endif
  26.  
  27. /*************   Dynamic Generation of Matrices   ************************/
  28.  
  29. dMatrix __vf  MD_matrix(  unsigned ht, unsigned len );
  30. dMatrix __vf  MD_matrix0( unsigned ht, unsigned len );
  31.     /*  notice that, in the memory model HUGE,
  32.         neither len nor ht may exceed 4095            */
  33.  
  34. /***************************************************************************
  35.  *  The following definitions ensure compatibility between dynamically     *
  36.  *  and statically allocated matrices. The definitions are somewhat        *
  37.  *  cumbersome, but the result for you is that you need not care about     *
  38.  *  the differences between the two types.                                 *
  39.  *  (Internally, the address of the first element of any matrix is needed; *
  40.  *  the expression "MA[0]" is evaluated in a different way for both types, *
  41.  *  but yields in either case the correct address to be passed to the      *
  42.  *  function you wish to call.)                                            *
  43.  *  Only in the rare case that you need to pass the address of one of      *
  44.  *  these functions as an argument to another function, you have to use    *
  45.  *  the actual run-time functions defined further below. Be careful with   *
  46.  *  this: future development of compilers may allow us to avoid this un-   *
  47.  *  handy scheme of macros. So future versions of MatrixLib may no longer  *
  48.  *  use these run-time names.                                              *
  49.  ***************************************************************************/
  50.  
  51.  
  52. /***  Addressing single elements of dynamically allocated matrices: ******
  53.      These two functions are for compatibility with Pascal
  54.      (where elements of dynamically allocated matrices are not directly
  55.      accessible), and for getting around the pointer arithmetics bug in
  56.      some versions of Borland C++.                                         */
  57.  
  58. #define MD_Pelement( MA, ht, len, m, n ) MDPelement( MA[0], ht, len, m, n )
  59.                      /* returns a pointer to MA[m][n]. */
  60. #define MD_element( MA, ht, len, m, n ) *MDPelement( MA[0], ht, len, m, n )
  61.                      /* dereferenced pointer */
  62.  
  63.  /****************  Initialization  ***************************************
  64.  
  65.     To initialize all elements of a matrix with the same value,
  66.     or to perform arithmetic operations on all elements simultaneously,
  67.     refer to the functions of VectorLib, declared in <VDstd.h>, <VDmath.h>.
  68.     In order to use the VectorLib functions, utilize the feature that
  69.     the whole matrix occupies one contiguous area in memory: pass the
  70.     address of the first row to the desired vector function, the size
  71.     of the "vector" being ht * len.
  72.     For example, initialize all elements of the matrix MA with 1.0
  73.     (this is *NOT* the identity matrix)  by calling
  74.         VD_equ1( MA[0], ht * len );
  75. */
  76.  
  77. #define MD_equ0( MA, ht, len )            VD_equ0( MA[0], ((ui)ht)*len )
  78. #define MD_equ1( MA, len )                MDequ1( MA[0], len )
  79.                        /* this is the identity matrix */
  80. #define MD_outerprod( MA, X, Y, ht, len ) MDouterprod( MA[0], X, Y, ht, len )
  81.                        /* sizX=ht, sizY=len */
  82. #define MD_Row_equC( MA, ht, len, iRow, C ) \
  83.                                         MDRow_equC( MA[0], ht, len, iRow, C )
  84. #define MD_Col_equC( MA, ht, len, iCol, C ) \
  85.                                         MDCol_equC( MA[0], ht, len, iCol, C )
  86. #define MD_Dia_equC( MA, len, C )       MDDia_equC( MA[0], len, C )
  87.  
  88. #define MD_Row_equV( MA, ht, len, iRow, X ) \
  89.                                         MDRow_equV( MA[0], ht, len, iRow, X )
  90. #define MD_Col_equV( MA, ht, len, iCol, X ) \
  91.                                         MDCol_equV( MA[0], ht, len, iCol, X )
  92. #define MD_Dia_equV( MA, len, X )       MDDia_equV( MA[0], len, X )
  93.  
  94. #define MD_equM( MB, MA, ht, len )  VD_equV( MB[0], MA[0], (ui)(ht)*(len) )
  95.  
  96. #define MD_UequL( MA, len ) MDUequL( MA[0], len )
  97. #define MD_LequU( MA, len ) MDLequU( MA[0], len )
  98.          /* copy lower-diagonal elements into upper-diagonal
  99.            (or vice versa) by index-reflection, so as to
  100.            get a symmetric matrix    */
  101.  
  102.             /* data-type conversions:  */
  103. #define M_DtoF( MF, MD, ht, len ) V_DtoF( MF[0], MD[0], ((ui)ht)*len )
  104. #define M_FtoD( MD, MF, ht, len ) V_FtoD( MD[0], MF[0], ((ui)ht)*len )
  105. #define M_EtoD( MD, ME, ht, len ) V_EtoD( MD[0], ME[0], ((ui)ht)*len )
  106. #define M_DtoE( ME, MD, ht, len ) V_DtoE( ME[0], MD[0], ((ui)ht)*len )
  107.  
  108.             /* suitable windows for MF_spectrum: */
  109. #define MD_Hanning( MA, ht, len )  MDHanning( MA[0], ht, len )
  110. #define MD_Parzen( MA, ht, len )   MDParzen( MA[0], ht, len )
  111. #define MD_Welch( MA, ht, len )    MDWelch( MA[0], ht, len )
  112.  
  113. /********  Extracting a submatrix and copying a submatrix back  *********/
  114.  
  115. #define MD_submatrix( MSub, subHt, subLen, \
  116.                       MSrce, srceHt, srceLen, \
  117.                       firstRowInCol, sampInCol, firstColInRow, sampInRow ) \
  118.                MDsubmatrix(  MSub[0], subHt, subLen, \
  119.                              MSrce[0], srceHt, srceLen, \
  120.                              firstRowInCol, sampInCol, firstColInRow, sampInRow )
  121.  
  122. #define MD_submatrix_equM( MDest, destHt, destLen, \
  123.                            firstRowInCol, sampInCol, firstColInRow, sampInRow, \
  124.                            MSrce, srceHt, srceLen ) \
  125.                MDsubmatrix_equM(  MDest[0], destHt, destLen, \
  126.                              firstRowInCol, sampInCol, firstColInRow, sampInRow, \
  127.                              MSrce[0], srceHt, srceLen )
  128.  
  129. /*****   Extracting a single row or a single column or the diagonal  ******
  130.  *       and storing it into a vector                                     */
  131.  
  132. #define MD_Row_extract( Y, MA, ht, len, iRow ) \
  133.                                      MDRow_extract( Y, MA[0], ht, len, iRow )
  134. #define MD_Col_extract( Y, MA, ht, len, iCol ) \
  135.                                      MDCol_extract( Y, MA[0], ht, len, iCol )
  136. #define MD_Dia_extract( Y, MA, len ) MDDia_extract( Y, MA[0], len )
  137.  
  138.  
  139. /*****************    Basic arithmetic operations *********************
  140.                       performed on one single row,
  141.                       or one single column of any matrix,
  142.                       or on the diagonal of a square matrix
  143.  
  144.     Note: In contrast to the analogous VectorLib functions, the operations
  145.     are performed in-place, i.e. the input matrix itself is changed  */
  146.  
  147. #define MD_Row_addC( MA, ht, len, iRow, C ) \
  148.                                      MDRow_addC( MA[0], ht, len, iRow, C )
  149. #define MD_Col_addC( MA, ht, len, iCol, C ) \
  150.                                      MDCol_addC( MA[0], ht, len, iCol, C )
  151. #define MD_Dia_addC( MA, len, C )    MDDia_addC( MA[0], len, C )
  152.  
  153. #define MD_Row_addV( MA, ht, len, iRow, X ) \
  154.                                      MDRow_addV( MA[0], ht, len, iRow, X )
  155. #define MD_Col_addV( MA, ht, len, iCol, X ) \
  156.                                      MDCol_addV( MA[0], ht, len, iCol, X )
  157. #define MD_Dia_addV( MA, len, X )    MDDia_addV( MA[0], len, X )
  158.  
  159. #define MD_Row_subC( MA, ht, len, iRow, C ) \
  160.                                      MDRow_addC( MA[0], ht, len, iRow, (-C) )
  161. #define MD_Col_subC( MA, ht, len, iCol, C ) \
  162.                                      MDCol_addC( MA[0], ht, len, iCol, (-C) )
  163. #define MD_Dia_subC( MA, len, C )    MDDia_addC( MA[0], len, (-C) )
  164.  
  165. #define MD_Row_subV( MA, ht, len, iRow, X ) \
  166.                                      MDRow_subV( MA[0], ht, len, iRow, X )
  167. #define MD_Col_subV( MA, ht, len, iCol, X ) \
  168.                                      MDCol_subV( MA[0], ht, len, iCol, X )
  169. #define MD_Dia_subV( MA, len, X )    MDDia_subV( MA[0], len, X )
  170.  
  171. #define MD_Row_subrC( MA, ht, len, iRow, C ) \
  172.                                      MDRow_subrC( MA[0], ht, len, iRow, C )
  173. #define MD_Col_subrC( MA, ht, len, iCol, C ) \
  174.                                      MDCol_subrC( MA[0], ht, len, iCol, C )
  175. #define MD_Dia_subrC( MA, len, C )   MDDia_subrC( MA[0], len, C )
  176.  
  177. #define MD_Row_subrV( MA, ht, len, iRow, X ) \
  178.                                      MDRow_subrV( MA[0], ht, len, iRow, X )
  179. #define MD_Col_subrV( MA, ht, len, iCol, X ) \
  180.                                      MDCol_subrV( MA[0], ht, len, iCol, X )
  181. #define MD_Dia_subrV( MA, len, X )   MDDia_subrV( MA[0], len, X )
  182.  
  183. #define MD_Row_mulC( MA, ht, len, iRow, C ) \
  184.                                      MDRow_mulC( MA[0], ht, len, iRow, C )
  185. #define MD_Col_mulC( MA, ht, len, iCol, C ) \
  186.                                      MDCol_mulC( MA[0], ht, len, iCol, C )
  187. #define MD_Dia_mulC( MA, len, C )    MDDia_mulC( MA[0], len, C )
  188.  
  189. #define MD_Row_mulV( MA, ht, len, iRow, X ) \
  190.                                      MDRow_mulV( MA[0], ht, len, iRow, X )
  191. #define MD_Col_mulV( MA, ht, len, iCol, X ) \
  192.                                      MDCol_mulV( MA[0], ht, len, iCol, X )
  193. #define MD_Dia_mulV( MA, len, X )    MDDia_mulV( MA[0], len, X )
  194.  
  195. #define MD_Row_divC( MA, ht, len, iRow, C ) \
  196.                                      MDRow_divC( MA[0], ht, len, iRow, C )
  197. #define MD_Col_divC( MA, ht, len, iCol, C ) \
  198.                                      MDCol_divC( MA[0], ht, len, iCol, C )
  199. #define MD_Dia_divC( MA, len, C )    MDDia_divC( MA[0], len, C )
  200.  
  201. #define MD_Row_divV( MA, ht, len, iRow, X ) \
  202.                                       MDRow_divV( MA[0], ht, len, iRow, X )
  203. #define MD_Col_divV( MA, ht, len, iCol, X ) \
  204.                                       MDCol_divV( MA[0], ht, len, iCol, X )
  205. #define MD_Dia_divV( MA, len, X )     MDDia_divV( MA[0], len, X )
  206.  
  207. #define MD_Row_divrC( MA, ht, len, iRow, C ) \
  208.                                       MDRow_divrC( MA[0], ht, len, iRow, C )
  209. #define MD_Col_divrC( MA, ht, len, iCol, C ) \
  210.                                       MDCol_divrC( MA[0], ht, len, iCol, C )
  211. #define MD_Dia_divrC( MA, len, C )    MDDia_divrC( MA[0], len, C )
  212.  
  213. #define MD_Row_divrV( MA, ht, len, iRow, X ) \
  214.                                       MDRow_divrV( MA[0], ht, len, iRow, X )
  215. #define MD_Col_divrV( MA, ht, len, iCol, X ) \
  216.                                       MDCol_divrV( MA[0], ht, len, iCol, X )
  217. #define MD_Dia_divrV( MA, len, X )    MDDia_divrV( MA[0], len, X )
  218.  
  219.  
  220. /******  One-dimensional vector operations **********************
  221.          performed along all rows or all columns simultaneously,
  222.          or along the diagonal of a square matrix                */
  223.  
  224. #define MD_Rows_max( Y, MA, ht, len )     MDRows_max( Y, MA[0], ht, len )
  225. #define MD_Cols_max( Y, MA, ht, len )     MDCols_max( Y, MA[0], ht, len )
  226. #define MD_Dia_max(  MA, len )            MDDia_max(  MA[0], len )
  227. #define MD_Rows_min( Y, MA, ht, len )     MDRows_min( Y, MA[0], ht, len )
  228. #define MD_Cols_min( Y, MA, ht, len )     MDCols_min( Y, MA[0], ht, len )
  229. #define MD_Dia_min(  MA, len )            MDDia_min(  MA[0], len )
  230.  
  231. #define MD_Rows_absmax( Y, MA, ht, len )  MDRows_absmax( Y, MA[0], ht, len )
  232. #define MD_Cols_absmax( Y, MA, ht, len )  MDCols_absmax( Y, MA[0], ht, len )
  233. #define MD_Dia_absmax(  MA, len )         MDDia_absmax(  MA[0], len )
  234. #define MD_Rows_absmin( Y, MA, ht, len )  MDRows_absmin( Y, MA[0], ht, len )
  235. #define MD_Cols_absmin( Y, MA, ht, len )  MDCols_absmin( Y, MA[0], ht, len )
  236. #define MD_Dia_absmin(  MA, len )         MDDia_absmin(  MA[0], len )
  237.  
  238. #define MD_Rows_sum( Y, MA, ht, len )     MDRows_sum( Y, MA[0], ht, len )
  239. #define MD_Cols_sum( Y, MA, ht, len )     MDCols_sum( Y, MA[0], ht, len )
  240. #define MD_Dia_sum(  MA, len )            MDDia_sum(  MA[0], len )
  241. #define MD_Rows_prod( Y, MA, ht, len )    MDRows_prod( Y, MA[0], ht, len )
  242. #define MD_Cols_prod( Y, MA, ht, len )    MDCols_prod( Y, MA[0], ht, len )
  243. #define MD_Dia_prod(  MA, len )           MDDia_prod(  MA[0], len )
  244.  
  245. #define MD_Rows_runsum( MA, ht, len )     MDRows_runsum( MA[0], ht, len )
  246. #define MD_Cols_runsum( MA, ht, len )     MDCols_runsum( MA[0], ht, len )
  247. #define MD_Rows_runprod( MA, ht, len )    MDRows_runprod( MA[0], ht, len )
  248. #define MD_Cols_runprod( MA, ht, len )    MDCols_runprod( MA[0], ht, len )
  249.  
  250. #define MD_Rows_rotate( MA, ht, len, pos) MDRows_rotate( MA[0], ht, len, pos )
  251. #define MD_Cols_rotate( MA, ht, len, pos) MDCols_rotate( MA[0], ht, len, pos )
  252. #define MD_Rows_reflect( MA, ht, len )    MDRows_reflect( MA[0], ht, len )
  253. #define MD_Cols_reflect( MA, ht, len )    MDCols_reflect( MA[0], ht, len )
  254.  
  255. /********  Operations involving two rows or two colums of one matrix  *****/
  256.  
  257. #define MD_Rows_exchange( MA, ht, len, row1, row2 ) \
  258.                              MDRows_exchange( MA[0], ht, len, row1, row2 )
  259. #define MD_Cols_exchange( MA, ht, len, col1, col2 ) \
  260.                              MDCols_exchange( MA[0], ht, len, col1, col2 )
  261.  
  262. #define MD_Rows_add( MA, ht, len, destRow, srceRow ) \
  263.                              MDRows_add( MA[0], ht, len, destRow, srceRow )
  264. #define MD_Cols_add( MA, ht, len, destCol, srceCol ) \
  265.                              MDCols_add( MA[0], ht, len, destCol, srceCol )
  266.                          /* dest[i] += source[i]  */
  267.  
  268. #define MD_Rows_sub( MA, ht, len, destRow, srceRow ) \
  269.                              MDRows_sub( MA[0], ht, len, destRow, srceRow )
  270. #define MD_Cols_sub( MA, ht, len, destCol, srceCol ) \
  271.                              MDCols_sub( MA[0], ht, len, destCol, srceCol )
  272.                          /* dest[i] -= source[i]  */
  273.  
  274. #define MD_Rows_Cadd( MA, ht, len, destRow, srceRow, C ) \
  275.                           MDRows_Cadd( MA[0], ht, len, destRow, srceRow, C )
  276. #define MD_Cols_Cadd( MA, ht, len, destCol, srceCol, C ) \
  277.                           MDCols_Cadd( MA[0], ht, len, destCol, srceCol, C )
  278.                          /* dest[i] += C * source[i]  */
  279.  
  280. #define MD_Rows_lincomb( MA, ht, len, destRow, destC, srceRow, srceC ) \
  281.             MDRows_lincomb( MA[0], ht, len, destRow, destC, srceRow, srceC )
  282. #define MD_Cols_lincomb( MA, ht, len, destCol, destC, srceCol, srceC ) \
  283.             MDCols_lincomb( MA[0], ht, len, destCol, destC, srceCol, srceC )
  284.   /*  linear combination: dest[i] = destC * dest[i] + sourceC * source[i] */
  285.  
  286.  
  287. /*************************  Transposing a matrix **********************/
  288.  
  289. #define MD_transpose( MTr, MA, htTr, lenTr ) \
  290.              MDtranspose( MTr[0], MA[0], htTr, lenTr )
  291.        /*  dimensions htTr, lenTr are those of the transposed matrix,
  292.            not of the original!                */
  293.  
  294.  
  295. /************************ Matrix Arithmetics *************************/
  296.  
  297. #define MD_addM( MC, MA, MB, htA, lenA ) \
  298.                  VD_addV( MC[0], MA[0], MB[0], ((ui)htA)*lenA )
  299. #define MD_subM( MC, MA, MB, htA, lenA ) \
  300.                  VD_subV( MC[0], MA[0], MB[0], ((ui)htA)*lenA )
  301. #define MD_mulC( MB, MA, htA, lenA, C ) \
  302.                  VD_mulC( MB[0], MA[0], ((ui)htA)*lenA, C )
  303. #define MD_divC( MB, MA, htA, lenA, C ) \
  304.                  VD_divC( MB[0], MA[0], ((ui)htA)*lenA, C )
  305. #define MDs_addM( MC, MA, MB, htA, lenA, C ) \
  306.                  VDs_addV( MC[0], MA[0], MB[0], ((ui)htA)*lenA, C )
  307. #define MDs_subM( MC, MA, MB, htA, lenA, C ) \
  308.                  VDs_subV( MC[0], MA[0], MB[0], ((ui)htA)*lenA, C )
  309. #define MD_lincomb( MC, MA, MB, htA, lenA, CA, CB ) \
  310.                  VD_lincomb( MC[0], MA[0], MB[0], ((ui)htA)*lenA, CA, CB )
  311.  
  312. #define MD_mulV( Y, MA, X, htA, lenA ) \
  313.                     MDmulV( Y, MA[0], X, htA, lenA )
  314.                           /*  Y = MA * X.  sizX=lenA, sizY=htA
  315.                               both X and Y are column-vectors    */
  316. #define VD_mulM( Y, X, MA, sizX, lenA ) \
  317.                     VDmulM( Y, X, MA[0], sizX, lenA )
  318.                          /*  Y = X * MA.  htA=sizX, sizY=lenA
  319.                              both X and Y are row-vectors.
  320.                              Mind the prefix: VD_ (not MD_)    */
  321. #define MD_mulM( MC, MA, MB, htA, lenA, lenB ) \
  322.                     MDmulM( MC[0], MA[0], MB[0], htA, lenA, lenB )
  323.                          /*  MC = MA * MB.  htB=lenA, htC=htA, lenC=lenB */
  324.  
  325.  
  326. /*************************  Linear Algebra    *****************************/
  327.  
  328.    /*  The standard treatment of linear systems is based
  329.        on LUD (matrix decomposition into Upper-triangular
  330.        and Lower-triangular components). The result of the
  331.        decomposition step is used for further operations.  */
  332.  
  333. #define MD_LUdecompose( MLU, Ind, MA, len ) \
  334.                     MDLUdecompose( MLU[0], Ind, MA[0], len )
  335.             /* returns "permut" = ±1 which is needed for MD_LUdet.
  336.                for singularities not cured by editing, permut is 0  */
  337. int     __vf MD_LUDresult( void );
  338.     /* returns 0, if MD_LUdecompose was successful;
  339.        returns 1, if MA was (nearly) singular in MD_LUdecompose.   */
  340. void    __vf MD_LUDsetEdit( double Thresh );
  341. double  __vf MD_LUDgetEdit( void );
  342.      /*  Editing threshold valid for MD_LUdecompose;
  343.          may be used to cure singularities           */
  344.  
  345. #define MD_LUsolve( X, MLU, B, Ind, len ) \
  346.                     MDLUsolve( X, MLU[0], B, Ind, len )
  347. #define MD_LUinv( MInv, MLU, Ind, len ) \
  348.                     MDLUinv( MInv[0], MLU[0], Ind, len )
  349. #define MD_LUdet( MLU, len, permut )  MDLUdet( MLU[0], len, permut )
  350. #define MD_LUimprove( X, B, MA, MLU, Ind, len ) \
  351.                       MDLUimprove( X, B, MA[0], MLU[0], Ind, len )
  352.  
  353.     /****  Special treatment of over- or under-determined
  354.            linear systems, i.e. of matrices with len != ht
  355.            and of singular matrices:
  356.            SVD (Singular Value Decomposition)       ****/
  357.  
  358. #define MD_SVdecompose( MU, MV, W, MA, htA, lenA ) \
  359.                     MDSVdecompose( MU[0], MV[0], W, MA[0], htA, lenA )
  360.             /*  sizB = htA,  sizX = sizW = htV = lenV = lenA  */
  361. #define MD_SVsolve( X, MU, MV, W, B, htU, lenU ) \
  362.                     MDSVsolve( X, MU[0], MV[0], W, B, htU, lenU )
  363.             /*  lenU = lenA,  htU = max( lenA, htA ) as fed into
  364.                 MD_SVdecompose   */
  365. #define MD_SVimprove( X, B, MA, MU, MV, W, htA, lenA ) \
  366.                     MDSVimprove( X, B, MA[0], MU[0], MV[0], W, htA, lenA )
  367. void    __vf MD_SVDsetEdit( double Thresh );
  368. double  __vf MD_SVDgetEdit( void );
  369.     /* Override of the standard values for editing threshholds
  370.        in MD_SVsolve. Calling MD_setEdit with Thresh=0.0 means
  371.        that you do the necessary editing of W yourself
  372.        before calling MD_SVsolve                           */
  373.  
  374.  /*****  "Easy-to-use" versions of the matrix functions
  375.           using LUD or SVD.
  376.           They allocate their own working space and rely
  377.           on your setting of the editing threshold. In
  378.           case of memory stress, you might better use the
  379.           two-step methods declared above.            ***/
  380. #define MD_solve( X, MA, B, len ) \
  381.                     MDsolve( X, MA[0], B, len )
  382. #define MD_inv( MInv, MA, len ) \
  383.                     MDinv( MInv[0], MA[0], len )
  384. #define MD_det(  MA, len ) \
  385.                     MDdet( MA[0], len )
  386. #define MD_solveBySVD( X, MA, B, ht, len ) \
  387.                     MDsolveBySVD( X, MA[0], B, ht, len )
  388.          /*  sizX = len,  sizB = ht  */
  389. #define MD_safeSolve( X, MA, B, len ) \
  390.                     MDsafeSolve( X, MA[0], B, len )
  391.         /* MD_safeSolve tries first LUD. If that fails, SVD is done.
  392.            X[i] will be 0.0 instead of INF for those i corresponding
  393.            to singularities. If even SVD fails, all X[i] are set to 0.0.
  394.            return value 0: success via LUD; 1: success via SVD; -1: failure */
  395.  
  396.       /*********  Eigenvalues and Eigenvectors  ********/
  397.       /*** only the most frequent case of symmetric real matrices
  398.            is covered here! *********/
  399.  
  400. #define MDsym_eigenvalues( EigV, EigM, MA, len, CalcEigenVec ) \
  401.                     MDs_eigenvalues( EigV, EigM[0], MA[0], len, CalcEigenVec )
  402.                         /*  Eigenvalues are returned in EigV,
  403.                             Eigenvectors are returned as the columns of EigM.
  404.                             CalcEigenVec = 0 means that only eigenvalues
  405.                             are needed; CalcEigenVec != 0 means that
  406.                             also eigenvectors are calculated.
  407.                             Even if eigenvectors are not desired, EigM is
  408.                             needed by the function as working-space. Then, on
  409.                             output, it will contain just rubbish.
  410.                             MA may be overwritten by EigM. */
  411.  
  412. /***************** Two-Dimensional Fourier-Transform Methods ************/
  413.  
  414. #if defined __cplusplus && defined _CMATH_CLASSDEFS
  415. } // following function cannot be extern "C"
  416. #endif
  417. dComplex  __vf   VD_getRspEdit( void );
  418. #if defined __cplusplus && defined _CMATH_CLASSDEFS
  419. extern "C" {
  420. #endif
  421. void      __vf   VD_setRspEdit( dComplex Trunc );
  422.            /* these are the same functions as used
  423.               in the one-dimensional case */
  424.  
  425. #define MDl_FFT( MY, MX, ht, len, dir )     MDlFFT( MY[0], MX[0], ht, len, dir )
  426. #define MDs_FFT( MY, MX, ht, len, dir )     MDsFFT( MY[0], MX[0], ht, len, dir )
  427. #define MDl_convolve( MY, MFlt, MX, MRsp, ht, len ) \
  428.                   MDlconvolve( MY[0], MFlt[0], MX[0], MRsp[0], ht, len )
  429. #define MDl_deconvolve( MY, MFlt, MX, MRsp, ht, len ) \
  430.                   MDldeconvolve( MY[0], MFlt[0], MX[0], MRsp[0], ht, len )
  431. #define MDl_filter( MY, MX, MFlt, ht, len ) MDlfilter( MY[0], MX[0], MFlt[0], ht, len )
  432. #define MDl_autocorr( MACorr, MX, ht, len ) MDlautocorr( MACorr[0], MX[0], ht, len )
  433. #define MDl_xcorr( MXCorr, MX, MY, ht, len) MDlxcorr( MXCorr[0], MX[0], MY[0], ht, len )
  434. #define MDl_spectrum( MSpec, htSpec, lenSpec, MX, htX, lenX, MWin ) \
  435.               MDlspectrum( MSpec[0], htSpec, lenSpec, MX[0], htX, lenX, MWin[0] )
  436.               /*   htSpec, lenSpec must be 2**n,
  437.                    MSpec must be a (htSpec+1)*(lenSpec+1) matrix!!
  438.                    htX >= n*htSpec,  lenX >= n*lenSpec,
  439.                    htWin = 2*htSpec, lenWin = 2*lenSpec */
  440.  
  441. #define MDs_convolve( MY, MFlt, MX, MRsp, ht, len ) \
  442.                   MDsconvolve( MY[0], MFlt[0], MX[0], MRsp[0], ht, len )
  443. #define MDs_deconvolve( MY, MFlt, MX, MRsp, ht, len ) \
  444.                   MDsdeconvolve( MY[0], MFlt[0], MX[0], MRsp[0], ht, len )
  445. #define MDs_filter( MY, MX, MFlt, ht, len ) MDsfilter( MY[0], MX[0], MFlt[0], ht, len )
  446. #define MDs_autocorr( MACorr, MX, ht, len ) MDsautocorr( MACorr[0], MX[0], ht, len )
  447. #define MDs_xcorr( MXCorr, MX, MY, ht, len) MDsxcorr( MXCorr[0], MX[0], MY[0], ht, len )
  448. #define MDs_spectrum( MSpec, htSpec, lenSpec, MX, htX, lenX, MWin ) \
  449.               MDsspectrum( MSpec[0], htSpec, lenSpec, MX[0], htX, lenX, MWin[0] )
  450.                /*  htSpec, lenSpec must be 2**n,
  451.                    MSpec has [htSpec+1][lenSpec+1] elements (!)
  452.                    htX >= n*htSpec,  lenX >= n*lenSpec,
  453.                    htWin = 2*htSpec, lenWin = 2*lenSpec     */
  454.  
  455. #if defined( __LARGE__ ) || defined( __COMPACT__ ) || defined ( __HUGE__ )
  456.    #define MD_FFT         MDl_FFT
  457.    #define MD_convolve    MDl_convolve
  458.    #define MD_deconvolve  MDl_deconvolve
  459.    #define MD_filter      MDl_filter
  460.    #define MD_autocorr    MDl_autocorr
  461.    #define MD_xcorr       MDl_xcorr
  462.    #define MD_spectrum    MDl_spectrum
  463. #else
  464.    #define MD_FFT         MDs_FFT
  465.    #define MD_convolve    MDs_convolve
  466.    #define MD_deconvolve  MDs_deconvolve
  467.    #define MD_filter      MDs_filter
  468.    #define MD_autocorr    MDs_autocorr
  469.    #define MD_xcorr       MDs_xcorr
  470.    #define MD_spectrum    MDs_spectrum
  471. #endif
  472.  
  473.  
  474. /************************** Data Fitting *********************************
  475.  
  476.   Notice that some of these functions have the prefix VD_, others MD_.
  477.   This depends on the form in which the data to be fitted are recorded:
  478.   vectors are fitted by the VD_ functions, matrices by the MD_ functions.
  479.   All of these functions employ matrix methods internally. The weighted
  480.   versions return covariances in a matrix "Covar". So they are all
  481.   contained in MatrixLib and declared here.
  482. */
  483.  
  484. void __vf VD_polyfit( dVector ParValues, unsigned deg, dVector X, dVector Y, ui sizex );
  485. #define   VD_polyfitwW( ParValues, Covar, deg, X, Y, InvVar, sizex ) \
  486.           VDpolyfitwW( ParValues, Covar[0], deg, X, Y, InvVar, sizex )
  487.        /* the size of A is deg+1 and Covar has [deg+1][deg+1] elements! */
  488.  
  489. void __vf VD_linfit( dVector ParValues, iVector ParStatus, unsigned npars,
  490.                      dVector X, dVector Y, ui sizex,
  491.                      void (*funcs)(dVector BasFuncs, double x, unsigned nfuncs) );
  492. #define   VD_linfitwW( ParValues, Covar, ParStatus, npars, X, Y, InvVar, sizex, funcs ) \
  493.           VDlinfitwW( ParValues, Covar[0], ParStatus, npars, X, Y, InvVar, sizex, funcs )
  494. #define   MD_linfit( ParValues, ParStatus, npars, X, Y, MZ, htZ, lenZ, funcs ) \
  495.           MDlinfit( ParValues, ParStatus, npars, X, Y, MZ[0], htZ, lenZ, funcs )
  496. #define   MD_linfitwW( ParValues, Covar, ParStatus, npars, X, Y, MZ, MInvVar, htZ, lenZ, funcs ) \
  497.           MDlinfitwW( ParValues, Covar[0], ParStatus, npars, X, Y, MZ[0], MInvVar[0], htZ, lenZ, funcs )
  498.  
  499. void  __vf VD_setLinfitNeglect( double Thresh );
  500.                    /* neglect A[i]=0, if significance smaller than Thresh */
  501. double __vf VD_getLinfitNeglect( void );
  502.  
  503. double __vf VD_nonlinfit( dVector ParValues, iVector ParStatus, unsigned npars,
  504.                          dVector X, dVector Y, ui sizex,
  505.                          void (*modelfunc)(dVector YModel, dVector XModel, ui size),
  506.                          void (*derivatives)(dVector dYdPari,dVector X, ui size, unsigned i) );
  507.             /* returns figure-of-merit of best A. If you don't know the partial
  508.                derivatives with respect to ParValues, call with derivatives=NULL */
  509. #define   VD_nonlinfitwW( ParValues, Covar, ParStatus, npars, X, Y, InvVar, sizex, modelfunc, deriv ) \
  510.           VDnonlinfitwW( ParValues, Covar[0], ParStatus, npars, X, Y, InvVar, sizex, modelfunc, deriv )
  511. #define   MD_nonlinfit( ParValues, ParStatus, npars, X, Y, MZ, htZ, lenZ, modelfunc, deriv ) \
  512.           MDnonlinfit( ParValues, ParStatus, npars, X, Y, MZ[0], htZ, lenZ, modelfunc, deriv )
  513. #define   MD_nonlinfitwW( ParValues, Covar, ParStatus, npars, X, Y, MZ, MInvVar, htZ, lenZ, modelfunc, deriv ) \
  514.           MDnonlinfitwW( ParValues, Covar[0], ParStatus, npars, X, Y, MZ[0], MInvVar[0], htZ, lenZ, modelfunc, deriv )
  515.  
  516.         /* If you know some partial derivatives, you may call these functions
  517.            for those parameters for which you do not know them:           */
  518. void   __vf VD_nonlinfit_autoDeriv( dVector dYdPari, dVector X, ui size, unsigned ipar );
  519. void   __vf VD_nonlinfitwW_autoDeriv( dVector dYdPari, dVector X, ui size, unsigned ipar );
  520. void   __vf MD_nonlinfitwW_autoDeriv( dMatrix dZdAi, unsigned htZ, unsigned lenZ, dVector X, dVector Y, unsigned ipar );
  521. void   __vf MD_nonlinfit_autoDeriv( dMatrix dZdAi, unsigned htZ, unsigned lenZ, dVector X, dVector Y, unsigned ipar );
  522.        /* The following functions allow to monitor the progress of
  523.           a nonlinear fitting procedure or to manually stop it:     */
  524. double   __vf VD_nonlinfit_getChi2( void );
  525. double   __vf VD_nonlinfitwW_getChi2( void );
  526. double   __vf MD_nonlinfit_getChi2( void );
  527. double   __vf MD_nonlinfitwW_getChi2( void );
  528. void     __vf VD_nonlinfit_getBestValues( dVector BestValues );
  529. void     __vf VD_nonlinfitwW_getBestValues( dVector BestValues );
  530. void     __vf MD_nonlinfit_getBestValues( dVector BestValues );
  531. void     __vf MD_nonlinfitwW_getBestValues( dVector BestValues );
  532. unsigned __vf VD_nonlinfit_getTestRun( void );
  533. unsigned __vf VD_nonlinfitwW_getTestRun( void );
  534. unsigned __vf MD_nonlinfit_getTestRun( void );
  535. unsigned __vf MD_nonlinfitwW_getTestRun( void );
  536. unsigned __vf VD_nonlinfit_getTestPar( void );
  537. unsigned __vf VD_nonlinfitwW_getTestPar( void );
  538. unsigned __vf MD_nonlinfit_getTestPar( void );
  539. unsigned __vf MD_nonlinfitwW_getTestPar( void );
  540. int      __vf VD_nonlinfit_getTestDir( void );
  541. int      __vf VD_nonlinfitwW_getTestDir( void );
  542. int      __vf MD_nonlinfit_getTestDir( void );
  543. int      __vf MD_nonlinfitwW_getTestDir( void );
  544. void     __vf VD_nonlinfit_stop( void );
  545. void     __vf VD_nonlinfitwW_stop( void );
  546. void     __vf MD_nonlinfit_stop( void );
  547. void     __vf MD_nonlinfitwW_stop( void );
  548.  
  549. #ifdef __BORLANDC__
  550.     #pragma option -a-  /* avoid insertion of dummy bytes */
  551. #else   /* MS Visual C++ */
  552.     #pragma pack(push,1)
  553. #endif    /*  Borland or Microsoft */
  554.  
  555. typedef struct VD_NONLINFITOPTIONS
  556. {
  557.       int        FigureOfMerit;  /*  0:least squares, 1:robust */
  558.              /* Convergence conditions: if the changes achieved
  559.                 in successive iterations are lower than any of the
  560.                 following values, this signals convergence. Set
  561.                 criteria to 0.0, if not applicable              */
  562.       double     AbsTolChi,  /* absolute change of chi */
  563.                  FracTolChi, /* fractional change of chi */
  564.                  AbsTolPar,  /* absolute change of all parameters */
  565.                  FracTolPar; /* fractional change of all parameters */
  566.       unsigned   HowOftenFulfill; /* how often fulfill the above conditions? */
  567.       unsigned   LevelOfMethod;  /* 1: Levenberg-Marquardt method,
  568.                                     2: Downhill Simplex (Nelder and Mead) method,
  569.                                     3: both methods alternating;
  570.                                        add 4 to this in order to try
  571.                                        breaking out of local minima;
  572.                                     0: no fit, calculate only chi2 (and Covar) */
  573.       unsigned   LevMarIterations; /* max.number of successful iterations of LevMar */
  574.       unsigned   LevMarStarts; /* number of starts per LevMar run */
  575.       double     LambdaStart,
  576.                  LambdaMin, LambdaMax,
  577.                  LambdaDiv, LambdaMul;    /* LevMar parameter lambda */
  578.       unsigned   DownhillIterations; /* max.number of successful iterations in Downhill */
  579.       double     DownhillReflection,  /* re-shaping of the simplex */
  580.                  DownhillContraction,
  581.                  DownhillExpansion;
  582.       unsigned   TotalStarts;  /* max. number of LevMar/Downhill pairs */
  583.       dVector    UpperLimits;  /* impose upper limits on parameters */
  584.       dVector    LowerLimits;  /* impose lower limits on parameters */
  585.       void       (*Restrictions)(void);  /* user-defined editing of parameters */
  586. }  VD_NONLINFITOPTIONS;
  587.  
  588. void __vf VD_setNonlinfitOptions( VD_NONLINFITOPTIONS *Options );
  589. void __vf VD_getNonlinfitOptions( VD_NONLINFITOPTIONS *Options );
  590.  
  591. typedef struct VD_EXPERIMENT
  592. {
  593.     dVector  X, Y, InvVar;
  594.     ui       size;
  595.     double   WeightOfExperiment;
  596.      /* InvVar and WeightOfExperiment are needed only for the
  597.         weighted variants of the multifit functions */
  598. } VD_EXPERIMENT;
  599.  
  600. typedef struct MD_EXPERIMENT
  601. {
  602.     dVector    X, Y;
  603.     dMatrix    MZ, MInvVar;  /* here no compatibility with static matrices! */
  604.     unsigned   htZ, lenZ;
  605.     double     WeightOfExperiment;
  606.      /* MInvVar and WeightOfExperiment are needed only for the
  607.         weighted variants of the multifit functions */
  608. } MD_EXPERIMENT;
  609.  
  610. #ifdef __BORLANDC__
  611.     #pragma option -a.   /* restore default data packing  */
  612. #else   /* MS Visual C++ */
  613.     #pragma pack(pop)
  614. #endif
  615.  
  616. void __vf VD_multiLinfit( dVector ParValues, iVector ParStatus, unsigned nParameters,
  617.                 VD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nExperiments,
  618.                 void (*funcs)(dVector BasFuncs, double x,
  619.                               unsigned nfuncs, unsigned iexperiment) );
  620. #define   VD_multiLinfitwW( ParValues, Covar, ParStatus, npars, ListOfEx, nexp, funcs ) \
  621.           VDmultiLinfitwW( ParValues, Covar[0], ParStatus, npars, ListOfEx, nexp, funcs )
  622. void __vf MD_multiLinfit( dVector ParValues, iVector ParStatus, unsigned nParameters,
  623.                 MD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nExperiments,
  624.                 void (*funcs)(dVector BasFuncs, double x, double y,
  625.                               unsigned nfuncs, unsigned iexperiment) );
  626. #define   MD_multiLinfitwW( ParValues, Covar, ParStatus, npars, ListOfEx, nexp, funcs ) \
  627.           MDmultiLinfitwW( ParValues, Covar[0], ParStatus, npars, ListOfEx, nexp, funcs )
  628. double __vf VD_multiNonlinfit( dVector ParValues,
  629.                 iVector ParStatus, unsigned npars,
  630.                 VD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nExperiments,
  631.                 void (*modelfunc)(dVector YModel, dVector XModel,
  632.                                   ui size, unsigned iexperiment),
  633.                 void (*derivatives)(dVector dYdPari,dVector X, ui size,
  634.                                    unsigned ipar, unsigned iexperiment) );
  635.             /* returns figure-of-merit of best ParValues. Chi2Detail contains
  636.                the chi2 contribution of each experiment.
  637.                If you don't know the partial derivatives, set derivatives=NULL */
  638. #define   VD_multiNonlinfitwW( ParValues, Covar, ParStatus, npars, ListOfEx, nexp, modfunc, deriv ) \
  639.           VDmultiNonlinfitwW( ParValues, Covar[0], ParStatus, npars, ListOfEx, nexp, modfunc, deriv )
  640. double __vf MD_multiNonlinfit( dVector ParValues, iVector ParStatus, unsigned npars,
  641.                 MD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nexperiments,
  642.                 void (*modelfunc)(dMatrix MZModel, unsigned htZ, unsigned lenZ,
  643.                                   dVector X, dVector Y, unsigned iexperiment),
  644.                 void (*derivatives)(dMatrix dZdAi, unsigned htZ, unsigned lenZ,
  645.                                     dVector X, dVector Y, unsigned ipar,
  646.                                     unsigned iexperiment) );
  647. #define   MD_multiNonlinfitwW( ParValues, Covar, ParStatus, npars, ListOfEx, nexp, modfunc, deriv ) \
  648.           MDmultiNonlinfitwW( ParValues, Covar[0], ParStatus, npars, ListOfEx, nexp, modfunc, deriv )
  649.  
  650. void  __vf VD_multiNonlinfit_autoDeriv( dVector dYdPari, dVector X, ui size,
  651.                                        unsigned iexperiment, unsigned ipar );
  652. void  __vf VD_multiNonlinfitwW_autoDeriv( dVector dYdPari, dVector X, ui size,
  653.                                        unsigned iexperiment, unsigned ipar );
  654. void  __vf MD_multiNonlinfit_autoDeriv( dMatrix dZdAi, unsigned htZ, unsigned lenZ,
  655.                                     dVector X, dVector Y,
  656.                                     unsigned ipar, unsigned iexperiment );
  657. void  __vf MD_multiNonlinfitwW_autoDeriv( dMatrix dZdAi, unsigned htZ, unsigned lenZ,
  658.                                     dVector X, dVector Y,
  659.                                     unsigned ipar, unsigned iexperiment );
  660. double   __vf VD_multiNonlinfit_getChi2( void );
  661. double   __vf VD_multiNonlinfitwW_getChi2( void );
  662. double   __vf MD_multiNonlinfit_getChi2( void );
  663. double   __vf MD_multiNonlinfitwW_getChi2( void );
  664. void     __vf VD_multiNonlinfit_getChi2Detail( dVector Chi2Detail );
  665. void     __vf VD_multiNonlinfitwW_getChi2Detail( dVector Chi2Detail );
  666. void     __vf MD_multiNonlinfit_getChi2Detail( dVector Chi2Detail );
  667. void     __vf MD_multiNonlinfitwW_getChi2Detail( dVector Chi2Detail );
  668. void     __vf VD_multiNonlinfit_getBestValues( dVector BestValues );
  669. void     __vf VD_multiNonlinfitwW_getBestValues( dVector BestValues );
  670. void     __vf MD_multiNonlinfit_getBestValues( dVector BestValues );
  671. void     __vf MD_multiNonlinfitwW_getBestValues( dVector BestValues );
  672. unsigned __vf VD_multiNonlinfit_getTestRun( void );
  673. unsigned __vf VD_multiNonlinfitwW_getTestRun( void );
  674. unsigned __vf MD_multiNonlinfit_getTestRun( void );
  675. unsigned __vf MD_multiNonlinfitwW_getTestRun( void );
  676. unsigned __vf VD_multiNonlinfit_getTestPar( void );
  677. unsigned __vf VD_multiNonlinfitwW_getTestPar( void );
  678. unsigned __vf MD_multiNonlinfit_getTestPar( void );
  679. unsigned __vf MD_multiNonlinfitwW_getTestPar( void );
  680. int      __vf VD_multiNonlinfit_getTestDir( void );
  681. int      __vf VD_multiNonlinfitwW_getTestDir( void );
  682. int      __vf MD_multiNonlinfit_getTestDir( void );
  683. int      __vf MD_multiNonlinfitwW_getTestDir( void );
  684. void     __vf VD_multiNonlinfit_stop( void );
  685. void     __vf VD_multiNonlinfitwW_stop( void );
  686. void     __vf MD_multiNonlinfit_stop( void );
  687. void     __vf MD_multiNonlinfitwW_stop( void );
  688.  
  689. /**************************  Input and Output  **************************/
  690.  
  691. #define MD_fprint( stream, MA, ht, len, linewidth ) \
  692.                     MDfprint( stream, MA[0], ht, len, linewidth )
  693. #define MD_print( MA, ht, len )  MDfprint( stdout, MA[0], ht, len, 80 )
  694. #define MD_cprint( MA, ht, len ) MDcprint( MA[0], ht, len )
  695.   /*  MD_print, MD_cprint usable only for DOS, EasyWin, and Win32 console applications! */
  696.  
  697. #define MD_write( str, MA, ht, len )   MDwrite( str, MA[0], ht, len )
  698. #define MD_read( MA, ht, len, str )    MDread( MA[0], ht, len, str )
  699. #define MD_setWriteFormat              VD_setWriteFormat
  700. #define MD_setWriteSeparate            VD_setNWriteSeparate
  701.                                    /* write and read in ascii format */
  702. #define MD_store( str, MA, ht, len ) \
  703.                            VD_store( str, MA[0], ((ui)(len))*(ht) );
  704. #define MD_recall( MA, ht, len, str) \
  705.                            VD_recall( MA[0], ((ui)(len))*(ht), str);
  706.                                   /* store and recall in binary format */
  707.  
  708.  
  709. /*************************************************************************
  710.  * Here are now the actual declarations of the functions used internally.*
  711.  * These declarations may change in future versions of MatrixLib!        *
  712.  * You should not care too much about them, except in the case you need  *
  713.  * the actual address of a run-time function (see above). Under all      *
  714.  * "normal" circumstances, use only the names defined above in the       *
  715.  * macro section!                                                        *
  716.  *************************************************************************/
  717.  
  718. double _VFAR * MDPelement( dPMatrix X, unsigned ht, unsigned len,
  719.                          unsigned m, unsigned n );
  720.                   /* pointer is normalized in memory model HUGE */
  721.  
  722. #define MDequ0( MA, ht, len )            VD_equ0( MA, ((ui)ht)*len )
  723. void   __vf  MDequ1( dPMatrix MA, unsigned len );  /* identity matrix */
  724. void   __vf  MDouterprod( dPMatrix MA, dVector X,  dVector Y,
  725.                           unsigned ht, unsigned len );
  726.  
  727. void   __vf  MDRow_equC( dPMatrix MA, unsigned ht, unsigned len,
  728.                          unsigned iRow, double C );
  729. void   __vf  MDCol_equC( dPMatrix MA, unsigned ht, unsigned len,
  730.                          unsigned iCol, double C );
  731. void   __vf  MDDia_equC( dPMatrix MA, unsigned len, double C );
  732.  
  733. void   __vf  MDRow_equV( dPMatrix MA, unsigned ht, unsigned len,
  734.                          unsigned iRow, dVector X );
  735. void   __vf  MDCol_equV( dPMatrix MA, unsigned ht, unsigned len,
  736.                          unsigned iCol, dVector X );
  737. void   __vf  MDDia_equV( dPMatrix MA, unsigned len, dVector X );
  738.  
  739. #define MDequM( MB, MA, ht, len )  VD_equV( MB, MA, (ui)(ht)*(len) )
  740. void   __vf  MDUequL( dPMatrix MA, unsigned len );
  741. void   __vf  MDLequU( dPMatrix MA, unsigned len );
  742.  
  743. void   __vf  MDHanning( dPMatrix MA, unsigned ht, unsigned len );
  744. void   __vf  MDParzen(  dPMatrix MA, unsigned ht, unsigned len );
  745. void   __vf  MDWelch(   dPMatrix MA, unsigned ht, unsigned len );
  746.  
  747. /********  Extracting a submatrix and copying a submatrix back  *********/
  748.  
  749. void  __vf  MDsubmatrix( dPMatrix MSub,
  750.                           unsigned subHt,  unsigned subLen,
  751.                           dPMatrix MSrce,
  752.                           unsigned srceHt,  unsigned srceLen,
  753.                           unsigned firstRowInCol,  unsigned sampInCol,
  754.                           unsigned firstColInRow,  unsigned sampInRow );
  755.  
  756. void  __vf  MDsubmatrix_equM( dPMatrix MDest,
  757.                                unsigned destHt,     unsigned destLen,
  758.                                unsigned firstRowInCol,  unsigned sampInCol,
  759.                                unsigned firstColInRow,  unsigned sampInRow,
  760.                                dPMatrix MSrce,
  761.                                unsigned srceHt,     unsigned srceLen );
  762.  
  763.  
  764. /*****   Extracting a single row or a single column or the diagonal  *****
  765.  *       and storing it into a vector                                    */
  766.  
  767. void __vf MDRow_extract( dVector Y, dPMatrix MA, unsigned ht, unsigned len,
  768.                            unsigned iRow );
  769. void __vf MDCol_extract( dVector Y, dPMatrix MA, unsigned ht, unsigned len,
  770.                            unsigned iCol );
  771. void __vf MDDia_extract( dVector Y, dPMatrix MA, unsigned len );
  772.  
  773.  
  774. /*****************    Basic arithmetic operations ***********************
  775.      performed on one single row,  or one single column of any matrix,
  776.      or on the diagonal of a square matrix                              */
  777.  
  778. void   __vf  MDRow_addC( dPMatrix MA, unsigned ht, unsigned len,
  779.                          unsigned iRow, double C );
  780. void   __vf  MDCol_addC( dPMatrix MA, unsigned ht, unsigned len,
  781.                          unsigned iCol, double C );
  782. void   __vf  MDDia_addC( dPMatrix MA, unsigned len, double C );
  783.  
  784. void   __vf  MDRow_subC( dPMatrix MA, unsigned ht, unsigned len,
  785.                          unsigned iRow, double C );
  786. void   __vf  MDCol_subC( dPMatrix MA, unsigned ht, unsigned len,
  787.                          unsigned iCol, double C );
  788. void   __vf  MDDia_subC( dPMatrix MA, unsigned len, double C );
  789.  
  790. void   __vf  MDRow_addV( dPMatrix MA, unsigned ht, unsigned len,
  791.                          unsigned iRow, dVector X );
  792. void   __vf  MDCol_addV( dPMatrix MA, unsigned ht, unsigned len,
  793.                          unsigned iCol, dVector X );
  794. void   __vf  MDDia_addV( dPMatrix MA, unsigned len, dVector X );
  795.  
  796. void   __vf  MDRow_subV( dPMatrix MA, unsigned ht, unsigned len,
  797.                          unsigned iRow, dVector X );
  798. void   __vf  MDCol_subV( dPMatrix MA, unsigned ht, unsigned len,
  799.                          unsigned iCol, dVector X );
  800. void   __vf  MDDia_subV( dPMatrix MA, unsigned len, dVector X );
  801.  
  802. void   __vf  MDRow_subrC( dPMatrix MA, unsigned ht, unsigned len,
  803.                          unsigned iRow, double C );
  804. void   __vf  MDCol_subrC( dPMatrix MA, unsigned ht, unsigned len,
  805.                          unsigned iCol, double C );
  806. void   __vf  MDDia_subrC( dPMatrix MA, unsigned len, double C );
  807.  
  808. void   __vf  MDRow_subrV( dPMatrix MA, unsigned ht, unsigned len,
  809.                          unsigned iRow, dVector X );
  810. void   __vf  MDCol_subrV( dPMatrix MA, unsigned ht, unsigned len,
  811.                          unsigned iCol, dVector X );
  812. void   __vf  MDDia_subrV( dPMatrix MA, unsigned len, dVector X );
  813.  
  814. void   __vf  MDRow_mulC( dPMatrix MA, unsigned ht, unsigned len,
  815.                          unsigned iRow, double C );
  816. void   __vf  MDCol_mulC( dPMatrix MA, unsigned ht, unsigned len,
  817.                          unsigned iCol, double C );
  818. void   __vf  MDDia_mulC( dPMatrix MA, unsigned len, double C );
  819.  
  820. void   __vf  MDRow_mulV( dPMatrix MA, unsigned ht, unsigned len,
  821.                          unsigned iRow, dVector X );
  822. void   __vf  MDCol_mulV( dPMatrix MA, unsigned ht, unsigned len,
  823.                          unsigned iCol, dVector X );
  824. void   __vf  MDDia_mulV( dPMatrix MA, unsigned len, dVector X );
  825.  
  826. void   __vf  MDRow_divC( dPMatrix MA, unsigned ht, unsigned len,
  827.                          unsigned iRow, double C );
  828. void   __vf  MDCol_divC( dPMatrix MA, unsigned ht, unsigned len,
  829.                          unsigned iCol, double C );
  830. void   __vf  MDDia_divC( dPMatrix MA, unsigned len, double C );
  831.  
  832. void   __vf  MDRow_divV( dPMatrix MA, unsigned ht, unsigned len,
  833.                          unsigned iRow, dVector X );
  834. void   __vf  MDCol_divV( dPMatrix MA, unsigned ht, unsigned len,
  835.                          unsigned iCol, dVector X );
  836. void   __vf  MDDia_divV( dPMatrix MA, unsigned len, dVector X );
  837.  
  838. void   __vf  MDRow_divrC( dPMatrix MA, unsigned ht, unsigned len,
  839.                          unsigned iRow, double C );
  840. void   __vf  MDCol_divrC( dPMatrix MA, unsigned ht, unsigned len,
  841.                          unsigned iCol, double C );
  842. void   __vf  MDDia_divrC( dPMatrix MA, unsigned len, double C );
  843.  
  844. void   __vf  MDRow_divrV( dPMatrix MA, unsigned ht, unsigned len,
  845.                          unsigned iRow, dVector X );
  846. void   __vf  MDCol_divrV( dPMatrix MA, unsigned ht, unsigned len,
  847.                          unsigned iCol, dVector X );
  848. void   __vf  MDDia_divrV( dPMatrix MA, unsigned len, dVector X );
  849.  
  850.  
  851. /******  One-dimensional vector operations **********************
  852.          performed along all rows or all columns simultaneously,
  853.          or along the diagonal of a square matrix                */
  854.  
  855. void   __vf  MDRows_max( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  856. void   __vf  MDCols_max( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  857. double __vf  MDDia_max(  dPMatrix MA, unsigned len );
  858. void   __vf  MDRows_min( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  859. void   __vf  MDCols_min( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  860. double __vf  MDDia_min(  dPMatrix MA, unsigned len );
  861.  
  862. void   __vf  MDRows_absmax( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  863. void   __vf  MDCols_absmax( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  864. double __vf  MDDia_absmax(  dPMatrix MA, unsigned len );
  865. void   __vf  MDRows_absmin( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  866. void   __vf  MDCols_absmin( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  867. double __vf  MDDia_absmin(  dPMatrix MA, unsigned len );
  868.  
  869. void   __vf  MDRows_sum( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  870. void   __vf  MDCols_sum( dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  871. double __vf  MDDia_sum(  dPMatrix MA, unsigned len );
  872. void   __vf  MDRows_prod(dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  873. void   __vf  MDCols_prod(dVector Y, dPMatrix MA, unsigned ht, unsigned len );
  874. double __vf  MDDia_prod( dPMatrix MA, unsigned len );
  875.  
  876. void  __vf  MDRows_runsum( dPMatrix MA, unsigned ht, unsigned len );
  877. void  __vf  MDCols_runsum( dPMatrix MA, unsigned ht, unsigned len );
  878. void  __vf  MDRows_runprod( dPMatrix MA, unsigned ht, unsigned len );
  879. void  __vf  MDCols_runprod( dPMatrix MA, unsigned ht, unsigned len );
  880.  
  881. void  __vf  MDRows_rotate( dPMatrix MA, unsigned ht, unsigned len, int pos );
  882. void  __vf  MDCols_rotate( dPMatrix MA, unsigned ht, unsigned len, int pos );
  883. void  __vf  MDRows_reflect( dPMatrix MA, unsigned ht, unsigned len );
  884. void  __vf  MDCols_reflect( dPMatrix MA, unsigned ht, unsigned len );
  885.  
  886. /********  Operations involving two rows or two colums of one matrix  *****/
  887.  
  888. void   __vf  MDRows_exchange( dPMatrix MA, unsigned ht, unsigned len,
  889.                               unsigned i1, unsigned i2 );
  890. void   __vf  MDCols_exchange( dPMatrix MA, unsigned ht, unsigned len,
  891.                               unsigned i1, unsigned i2 );
  892.  
  893. void   __vf  MDRows_add( dPMatrix MA, unsigned ht, unsigned len,
  894.                           unsigned destRow, unsigned sourceRow );
  895. void   __vf  MDCols_add( dPMatrix MA, unsigned ht, unsigned len,
  896.                           unsigned destCol, unsigned sourceCol );
  897.  
  898. void   __vf  MDRows_sub( dPMatrix MA, unsigned ht, unsigned len,
  899.                           unsigned destRow, unsigned sourceRow );
  900. void   __vf  MDCols_sub( dPMatrix MA, unsigned ht, unsigned len,
  901.                           unsigned destCol, unsigned sourceCol );
  902.  
  903. void   __vf  MDRows_Cadd( dPMatrix MA, unsigned ht, unsigned len,
  904.                            unsigned destRow, unsigned sourceRow, double C );
  905. void   __vf  MDCols_Cadd( dPMatrix MA, unsigned ht, unsigned len,
  906.                            unsigned destCol, unsigned sourceCol, double C );
  907.  
  908. void   __vf  MDRows_lincomb( dPMatrix MA, unsigned ht, unsigned len,
  909.                               unsigned destRow,  double  destC,
  910.                               unsigned srceRow,  double  srceC );
  911. void   __vf  MDCols_lincomb( dPMatrix MA, unsigned ht, unsigned len,
  912.                               unsigned destCol,  double  destC,
  913.                               unsigned srceCol,  double  srceC );
  914.  
  915.  
  916. /*************************  Transposing a matrix **********************/
  917.  
  918. void  __vf  MDtranspose( dPMatrix MTr, dPMatrix MA,
  919.                           unsigned htTr, unsigned lenTr );
  920.  
  921.  
  922. /************************ Matrix Arithmetics *************************/
  923.  
  924. #define MDaddM( MC, MA, MB, htA, lenA ) \
  925.                  VD_addV( MC, MA, MB, ((ui)htA)*lenA )
  926. #define MDsubM( MC, MA, MB, htA, lenA ) \
  927.                  VD_subV( MC, MA, MB, ((ui)htA)*lenA )
  928. #define MDmulC( MB, MA, htA, lenA, C ) \
  929.                  VD_mulC( MB, MA, ((ui)htA)*lenA, C )
  930. #define MDdivC( MB, MA, htA, lenA, C ) \
  931.                  VD_divC( MB, MA, ((ui)htA)*lenA, C )
  932. #define MDsaddM( MC, MA, MB, htA, lenA, C ) \
  933.                  VDs_addV( MC, MA, MB, ((ui)htA)*lenA, C )
  934. #define MDssubM( MC, MA, MB, htA, lenA, C ) \
  935.                  VDs_subV( MC, MA, MB, ((ui)htA)*lenA, C )
  936. #define MDlincomb( MC, MA, MB, htA, lenA, CA, CB ) \
  937.                  VD_lincomb( MC, MA, MB, ((ui)htA)*lenA, CA, CB )
  938. void  __vf  MDmulV( dVector Y, dPMatrix MA, dVector X,
  939.                     unsigned htMA, unsigned lenA );
  940. void  __vf  VDmulM( dVector Y, dVector X, dPMatrix MA,
  941.                     unsigned sizX, unsigned lenA );
  942. void  __vf  MDmulM( dPMatrix MC, dPMatrix MA, dPMatrix MB,
  943.                       unsigned htMA, unsigned lenMA, unsigned lenB );
  944.  
  945.  
  946. /*************************  Linear Algebra    *****************************/
  947.  
  948. int    __vf  MDLUdecompose( dPMatrix MLU,  uiVector Ind, dPMatrix MA,
  949.                             unsigned len );
  950.  
  951. void   __vf  MDLUsolve( dVector X, dPMatrix MLU, dVector B, uiVector Ind,
  952.                         unsigned len );
  953. void   __vf  MDLUinv( dPMatrix MInv, dPMatrix MLU, uiVector Ind,
  954.                       unsigned len );
  955. double __vf  MDLUdet( dPMatrix MLU, unsigned len, int permut );
  956. void   __vf  MDLUimprove( dVector X, dVector B, dPMatrix MA, dPMatrix MLU,
  957.                           uiVector Ind, unsigned len );
  958.  
  959. int    __vf  MDSVdecompose( dPMatrix MU, dPMatrix MV, dVector W, dPMatrix MA,
  960.                            unsigned htMA, unsigned lenA );
  961. void   __vf  MDSVsolve( dVector X, dPMatrix MU, dPMatrix MV, dVector W,
  962.                        dVector B, unsigned htU, unsigned lenU );
  963. void   __vf  MDSVimprove(  dVector X, dVector B, dPMatrix MA,
  964.                        dPMatrix MU, dPMatrix MV, dVector W,
  965.                        unsigned htMA, unsigned lenA );
  966.  
  967.           /*  functions using LUD or SVD     */
  968. int    __vf  MDsolve( dVector X, dPMatrix MA, dVector B, unsigned len );
  969.                   /* ret.value != 0 signals error */
  970. int    __vf  MDinv( dPMatrix MInv, dPMatrix MA, unsigned len );
  971.                  /* ret.value != 0 signals error */
  972. double __vf  MDdet( dPMatrix MA, unsigned len );
  973.  
  974. int    __vf  MDsolveBySVD( dVector X, dPMatrix MA, dVector B,
  975.                            unsigned ht, unsigned len );
  976.               /*  sizX = lenMA,  sizB = htA.  ret.value != 0 signals failure */
  977. int    __vf  MDsafeSolve( dVector X, dPMatrix MA, dVector B, unsigned len );
  978.               /* ret.value 0: success via LUD; 1: success via SVD; -1: error */
  979.  
  980.        /*********  Eigenvalues and Eigenvectors  ********/
  981.  
  982. void __vf MDs_eigenvalues( dVector EigV, dPMatrix EigM, dPMatrix MA, unsigned len,
  983.                         int CalcEigenVec );
  984.  
  985. /*************  Two-Dimensional Fourier-Transform Methods *****************/
  986.  
  987. void  __vf   MDlFFT( dPMatrix MY, dPMatrix MX,
  988.                      unsigned ht, unsigned len, int dir );
  989. void  __vf   MDlconvolve( dPMatrix MY, dPMatrix MFlt, dPMatrix MX,
  990.                           dPMatrix MRsp, unsigned ht, unsigned len );
  991. void  __vf   MDldeconvolve( dPMatrix MY, dPMatrix MFlt, dPMatrix MX,
  992.                             dPMatrix MRsp, unsigned ht, unsigned len );
  993. void  __vf   MDlfilter( dPMatrix MY, dPMatrix MX, dPMatrix MFlt,
  994.                         unsigned ht, unsigned len );
  995. void  __vf   MDlautocorr( dPMatrix MACorr, dPMatrix MX,
  996.                           unsigned ht, unsigned len );
  997. void  __vf   MDlxcorr( dPMatrix MXCorr, dPMatrix MX, dPMatrix MY,
  998.                        unsigned ht, unsigned len );
  999. void  __vf   MDlspectrum( dPMatrix MSpec, unsigned htSpec, unsigned lenSpec,
  1000.                           dPMatrix MX, unsigned htX, unsigned lenX,
  1001.                           dPMatrix MWin );
  1002.  
  1003. void  __vf   MDsFFT( dPMatrix MY, dPMatrix MX,
  1004.                      unsigned ht, unsigned len, int dir );
  1005. void  __vf   MDsconvolve( dPMatrix MY, dPMatrix MFlt, dPMatrix MX,
  1006.                           dPMatrix MRsp, unsigned ht, unsigned len );
  1007. void  __vf   MDsdeconvolve( dPMatrix MY, dPMatrix MFlt, dPMatrix MX,
  1008.                             dPMatrix MRsp, unsigned ht, unsigned len );
  1009. void  __vf   MDsfilter( dPMatrix MY, dPMatrix MX, dPMatrix MFlt,
  1010.                         unsigned ht, unsigned len );
  1011. void  __vf   MDsautocorr( dPMatrix MACorr, dPMatrix MX,
  1012.                           unsigned ht, unsigned len );
  1013. void  __vf   MDsxcorr( dPMatrix MXCorr, dPMatrix MX, dPMatrix MY,
  1014.                        unsigned ht, unsigned len );
  1015. void  __vf   MDsspectrum( dPMatrix MSpec, unsigned htSpec, unsigned lenSpec,
  1016.                           dPMatrix MX, unsigned htX, unsigned lenX,
  1017.                           dPMatrix MWin );
  1018.  
  1019.       /***************  Data Fitting    ******************/
  1020.  
  1021. void __vf VDpolyfitwW( dVector ParValues, dPMatrix Covar, unsigned deg,
  1022.                         dVector X, dVector Y, dVector InvVar, ui sizex );
  1023. void __vf VDlinfitwW( dVector ParValues, dPMatrix Covar, iVector ParStatus, unsigned npars,
  1024.                     dVector X, dVector Y, dVector InvVar, ui sizex,
  1025.                     void (*funcs)(dVector BasFuncs, double x, unsigned nfuncs) );
  1026. double __vf VDnonlinfitwW( dVector ParValues, dPMatrix Covar, iVector ParStatus, unsigned npars,
  1027.                     dVector X, dVector Y, dVector InvVar, ui sizex,
  1028.                     void (*modelfunc)(dVector YModel, dVector X, ui size),
  1029.                     void (*derivatives)(dVector dYdPari, dVector X, ui size, unsigned i) );
  1030. void __vf MDlinfit( dVector ParValues, iVector ParStatus, unsigned npars,
  1031.                     dVector X, dVector Y, dPMatrix MZ, unsigned htZ, unsigned lenZ,
  1032.                     void (*funcs)(dVector BasFuncs, double x, double y, unsigned nfuncs) );
  1033. double __vf MDnonlinfit( dVector ParValues, iVector ParStatus, unsigned npars,
  1034.                     dVector X, dVector Y, dPMatrix MZ, unsigned htZ, unsigned lenZ,
  1035.                     void (*modelfunc)(dMatrix MZModel, unsigned htZ, unsigned lenZ, dVector X, dVector Y ),
  1036.                     void (*derivatives)(dMatrix dZdAi, unsigned htZ, unsigned lenZ, dVector X, dVector Y, unsigned i) );
  1037. void __vf MDlinfitwW( dVector ParValues, dPMatrix Covar, iVector ParStatus, unsigned npars,
  1038.                       dVector X, dVector Y, dPMatrix MZ, dPMatrix MInvVar, unsigned htZ, unsigned lenZ,
  1039.                       void (*funcs)(dVector BasFuncs, double x, double y, unsigned nfuncs) );
  1040. double __vf MDnonlinfitwW( dVector ParValues, dPMatrix Covar, iVector ParStatus, unsigned npars,
  1041.                     dVector X, dVector Y, dPMatrix MZ, dPMatrix MInvVar, unsigned htZ, unsigned lenZ,
  1042.                     void (*modelfunc)(dMatrix MZModel, unsigned htZ, unsigned lenZ, dVector X, dVector Y ),
  1043.                     void (*derivatives)(dMatrix dZdAi, unsigned htZ, unsigned lenZ, dVector X, dVector Y, unsigned i) );
  1044.  
  1045. void __vf VDmultiLinfitwW( dVector ParValues, dPMatrix Covar, iVector ParStatus, unsigned npars,
  1046.                 VD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nexperiments,
  1047.                 void (*funcs)(dVector BasFuncs, double x,
  1048.                               unsigned nfuncs, unsigned nexperiment) );
  1049. void __vf MDmultiLinfitwW( dVector ParValues, dPMatrix Covar,
  1050.                 iVector ParStatus, unsigned npars,
  1051.                 MD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nexperiments,
  1052.                 void (*funcs)(dVector BasFuncs, double x, double y,
  1053.                               unsigned nfuncs, unsigned nexperiment) );
  1054. double __vf VDmultiNonlinfitwW( dVector ParValues, dPMatrix Covar,
  1055.                 iVector ParStatus, unsigned npars,
  1056.                 VD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nexperiments,
  1057.                 void (*modelfunc)(dVector YModel, dVector X, ui size,
  1058.                                   unsigned iexperiment),
  1059.                 void (*derivatives)(dVector dYdPari, dVector X, ui size,
  1060.                                   unsigned ipar, unsigned iexperiment) );
  1061. double __vf MDmultiNonlinfitwW( dVector ParValues, dPMatrix Covar,
  1062.                 iVector ParStatus, unsigned npars,
  1063.                 MD_EXPERIMENT _VFAR *ListOfExperiments, unsigned nexperiments,
  1064.                 void (*modelfunc)(dMatrix MZModel, unsigned htZ, unsigned lenZ,
  1065.                                   dVector X, dVector Y, unsigned iexperiment ),
  1066.                 void (*derivatives)(dMatrix dZdAi, unsigned htZ, unsigned lenZ,
  1067.                                     dVector X, dVector Y,
  1068.                                     unsigned ipar, unsigned iexperiment) );
  1069.  
  1070.       /*************  Input and Output  ****************/
  1071.  
  1072. void __vf MDfprint( FILE _VFAR *stream, dPMatrix MA, unsigned ht,
  1073.                      unsigned len, unsigned linewidth );
  1074. void __vf MDcprint( dPMatrix MA, unsigned ht, unsigned len );
  1075. void  __vf    MDwrite( FILE _VFAR *stream, dPMatrix X, unsigned ht, unsigned len  );
  1076. void  __vf    MDread( dPMatrix X, unsigned ht, unsigned len, FILE _VFAR *stream );
  1077. #define MDstore( str, MA, ht, len ) \
  1078.                            VD_store( str, MA, ((ui)(len))*(ht) );
  1079. #define MDrecall( MA, ht, len, str) \
  1080.                            VD_recall( MA, ((ui)(len))*(ht), str);
  1081.  
  1082. #ifdef __cplusplus
  1083. }
  1084. #endif
  1085. #endif /* __MDSTD_H */
  1086.